前向投影的数值计算

要求

1. 编程实现扇形束或平行束CT前向投影数值计算（射线驱动、像素驱动、距离驱动，图像旋转+像素累加，或者其他投影方式），获取Shepp Logan仿体的前向投影, 显示对应的sinogram图，并将结果与上次作业中解析解获取的前向投影结果相对比，分析二者差异

2. 修改数值计算仿体图像的大小（比如256*256，512×512,1024×1024等），分析不同情形下计算的投影计算的差异

思路

​ 前向投影的实现方式有很多，例如上面提到的射线驱动、像素驱动、距离驱动、图像旋转累加等方式，其中图像旋转加累加的方式计算Sheep Logan仿体的前向投影最为简便。

​ 除了使用图像旋转累加的方法计算sinogram图，本文还实现了射线驱动的方法计算前向投影。

图像旋转+像素累加

可以拆分成两个步骤：

1. 使用双线性插值的方法旋转图像；
2. 沿探测器垂线方向累加像素值；

由sinogram图可以看到，先旋转图像再累加的方式能够较好的得到前向投影图

clc,clear
% 生成 Shepp Logan仿体图像
P = phantom(128);  % 生成128*128像素值的图像
K = imrotate(P,45,'bilinear','crop');
figure
subplot(1,2,1)
imshow(P)
title('Original image')
subplot(1,2,2)
imshow(K)
title('Rotate theta 45°')


% 参数设定
NumDetector = 128;      % 探测器数量
NumAngle = 360;         % 旋转角度

Rdata = zeros(NumDetector,NumAngle);    
for i=1:NumAngle
    K = imrotate(P,i,'bilinear','crop');
    Rdata(:,i) = sum(K,1);
end

figure
imshow(Rdata,[])
colorbar()
colormap(gray)
title("Shepp-Logan数值法投影图")

射线驱动

射线驱动方式可以分为以下几步：

1. 旋转探测器并计算每个探测晶体中心的垂线
2. 计算该垂线穿过图像中的像素，并累加像素值

由投影结果可以看到，在某些特定角度如90°和270°，前向投影结果较差，原因是在这些位置k = 0，函数表达式y = b。

clc,clear
% 生成 Shepp Logan仿体图像
P = phantom(128); 
figure
imshow(P)
title('Shepp-Logan')

% 参数设定
NumDetector = 150;      % 探测器数量
NumAngle = 360;         % 旋转角度


Rdata = zeros(NumDetector,NumAngle);    
for i=1:NumAngle
    for j =1:NumDetector
        [k,b] = CalSlope(size(P),i,j,NumDetector);
        Rdata(j,i) = CalDis(k,b,P);
    end 
end

figure
imshow(Rdata,[])
colorbar()
colormap(gray)
title("射线驱动")

function [k,b] = CalSlope(ImgSize, angle, num, NumDetector)
% 本函数用于计算探测器在特定偏转角度，穿过特定探测晶体中心的直线方程
% Imgsize 图像大小
% angle 偏转角度
% num 晶体排序
% NumberDetector 晶体总数

DisAbord = 100;
%角度弧度制转换
anglePi = pi*angle / 180;
%计算k
%k1 = tan(anglePi);
if angle == 0
    k = 1;
elseif angle == 90
    k = 0;
else
    k = -1/tan(anglePi);
end
%计算最中心的传感器的坐标
center_x = ImgSize(2)/2+sin(anglePi)*DisAbord;
center_y = ImgSize(1)/2-cos(anglePi)*DisAbord;
%算出第num个传感器距离中心的距离
DisCenter = (num-(NumDetector+1)/2);
%计算第num个传感器的坐标
x = center_x + DisCenter*cos(anglePi);
y = center_y + DisCenter*sin(anglePi);
%根据y=kx+b计算b
b = y - k*x;

end

%% function
function value = CalDis(k,b,Img)
value = 0;
for i = 1:size(Img,2)
    for j = 1:size(Img,1)
        dis = abs(k*(i-0.5)-(j-0.5)+b)/sqrt(k^2+1);
        if k == 0||k == 1
            judge = 1/2;
        else
            judge = 1/2*abs(k/sqrt(k^2+1));
        end
        if dis < judge
            value = value + Img(j,i);
        end
    end
end
end

结果

投影数值法与解析法比较

​ 从图中可以看出解析法的结果更为平滑，过度较为自然。在数值法投影图中能够较为明显的看到投影图的缝隙，在细节上有所缺失。

不同图像尺寸的比较

在数值法求解投影图中，使用phantom(N)创建特定尺寸大小的灰度图

P = phantom(128);  % 生成128*128像素值的图像

​ 在图像中可以看出，256*256图像的sinogram图较128*128投影图有比较明显的改进，在细节上更加细腻，并且图像colorbar更长，所表达的灰度梯度更多。